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The  purpose  of  this  document  is  to  summarize  the  work  accomplished  under  work  unit  Q02Z.  The 
sole  source  of  funding  for  this  program  was  through  AFOSR  and  was  supported  by  Dr.  Jason 
Marshall  under  contract  14RQ05COR  and  Dr.  Fariba  Fahroo  under  contract  12RZ06COR. 

The  overarching  goals  of  this  project  was  to  develop  a  new  numerical  approach  for  multiscale 
plasma  simulations  that  spans  the  range  of  far- from  equilibrium  to  local  thermodynamic  equilibrium 
conditions.  This  is  a  very  broad  and  ambitious  goal  and,  in  order  to  bring  focus  to  this  project, 
the  plasma  parameters  were  limited  to  those  found  in  non-relativistic  laser  plasma  interactions. 
The  use  of  a  detailed  Collisional  Radiative  operator  was  necessary  to  accurately  account  for  all 
processes  of  the  collisional  cascade  during  the  relaxation  of  a  hot  plasma.  To  this  end,  the  focus 
was  in  the  development  of 

(a)  a  detailed  time  accurate  Collisional  Radiative  model,  (b)  a  one- dimensional  Multi-Fluid  for¬ 
mulation  that  includes  inelastic  collisions,  (c)  a  conservative  BGK  collision  operator,  and  (d)  hy¬ 
bridization  techniques  that  blends  low  and  high-fidelity  algorithms. 

The  work  performed  over  the  last  three  years  has  done  much  to  advance  the  state  of  the  art 
in  complexity  reduction  of  Collisional  Radiative  (CR)  modeling  as  well  as  the  effects  of  inelastic 
collisions  on  the  Multi-Fluid  description  of  plasmas.  This  work  has  been  recognized  in  two  workshop 
presentations,  two  conferences  papers  and  five  journal  articles;  four  published  [5,  6,  3,  7],  and  one 
awaiting  review.  The  remainder  of  this  report  will  summarize  the  accomplishments  and  techniques 
developed  during  this  investigation. 


1  Elastic  Collisions  and  CR  Models 

Argon  Collisional  Radiative  Modeling 

The  bulk  of  the  work  performed  over  the  last  three  years  focused  on  the  extension  of  the  CR 
model  originally  generated  under  a  prior  LRIR  Task.  In  the  previous  LRIR  task,  a  CR  model 
with  complexity  reduction  was  derived  and  verified  for  atomic  hydrogen.  This  investigation  focuses 
on  the  extension  of  the  previously  derived  algorithms  to  atoms  with  multiple  ionization  levels. 
This  work  was  split  between  two  successive  goals.  The  first  was  to  extend  the  current  model 
to  include  multiple  ionization  levels  without  complexity  reduction;  this  model  was  then  validated 
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1  ELASTIC  COLLISIONS  AND  CR  MODELS 


Figure  1:  Pressure  Vs.  Breakdown  intensity  comparison  with  experimental  results. 


against  experimental  results  as  well  as  with  code-to-code  comparison.  The  second  goal  was  to 
introduce  the  complexity  reduction  utilizing  both  uniform  and  Boltzmann  grouping  techniques. 
The  Boltzmann  and  uniform  grouping  techniques  were  then  compared  against  one  another.  The 
algorithms  and  corresponding  code  are  generic  and  valid  for  all  single  atomic  species  but  argon 
was  used  for  all  of  the  test  cases.  The  choice  of  argon  was  due  to  the  availability  of  a  complete 
set  of  collisional  cross-section  and  radiative  data,  made  available  through  the  Los  Alamos  National 
Labs  (LANL)  database[2,  4].  Argon  presence  in  a  wide  range  of  interests  (e.g.,  Field  Reverse 
Configuration  Thruster,  LPI  experiments)  was  also  a  contributing  factor  in  its  selection. 

For  the  experimental  comparison  the  code  was  compared  against  experiments  from  Sircar  [8],  shown 
in  Figure  1.  In  this  experiment  the  laser  power  required  for  breakdown  to  occur  was  measured  un¬ 
der  various  background  pressures.  The  main  difficulty  encountered  during  simulation  was  to  define 
when  breakdown  occurred  during  the  simulation  and  correlating  the  results  to  the  experimentally 
determined  breakdown  values.  Historically,  a  peak  ionization  fraction  of  0.1%  is  used.  Adopting 
the  0.1%  convention  resulted  in  good  agreement  for  the  532  nm  laser  case.  While  a  similar  trend 
was  observed  for  the  1064  nm  case,  the  numerical  result  deviates  for  higher  pressure  scenarios. 
This  deviation  was  most  likely  due  to  ionization  potential  depression  or  continuum  lowering  which 
was  not  accounted  for  in  the  code.  This  simulation  also  included  the  first  two  ionization  levels 
of  argon  for  a  total  of  633  excited  states,  101  for  Ar,  194  for  Ar+  and  338  Ar++.  To  accurately 
capture  the  various  CR  processes  required  modeling  several  elementary  processes  including  electron 
impact  excitation/de-excitation,  electron  impact  ionization/recombination,  elastic  collisions,  radi¬ 
ation  losses  due  to  bound-bound  transitions  and  bound-free  processes,  and  laser  based  processes 
including  multi-photon  ionization  and  inverse  Bremsstrahlung.  One  key  physical  process  not  in¬ 
cluded  is  radiation  transfer;  instead  an  optical  escape  factor  is  used  to  model  the  loss  of  energy  due 
to  radiation. 

The  second  benchmark  case  compares  our  code  against  FLYCHK,  a  CR  code  created  at  Lawrence 
Livermore  National  Labs[l].  The  ionization  level  vs  the  electron  temperature  was  used  to  compare 
the  two  codes  as  shown  in  Figure  2.  The  electron  temperature  varied  from  1  to  1000  eV  and  three 
different  number  densities,  10,  1020,  and  1023  cm-3,  were  used.  Overall,  the  two  codes  produced 
similar  results.  The  largest  discrepancies  occurred  with  the  high  density  case  at  low  temperatures 
and  the  low  density  case  at  high  temperature.  For  the  high  density  low  temperature  case,  the 
difference  is  likely  due  to  corrections  added  to  FLYCHK  that  account  for  high  density  effects  such 
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1  ELASTIC  COLLISIONS  AND  CR  MODELS 


Figure  2:  Comparison  of  equilibrium  average  ionization  between  FLYCHK  (dashed)  and  AFRL  (solid)  argon  CR 
solvers. 


as  continuum  lowering.  For  the  high  temperature  low  density  case,  radiative  recombination  rates 
dominate  three-body  collisional  deexcitation  which  establishes  equilibrium,  which  then  increases 
the  sensitivity  of  the  cross  section  data.  Because  the  two  codes  use  different  sets  of  atomic  data, 
discrepancies  in  this  regime  are  not  unexpected.  It  is  also  important  to  note  that  this  is  a  code-to- 
code  comparison  and  the  validation  of  these  simulations  is  an  area  of  active  research  in  the  non-local 
thermodynamic  equilibrium  physics  community. 

The  change  from  atomic  hydrogen  to  argon  required  updating  the  code  base  to  allow  for  non- 
analytical  cross  sections.  During  the  initial  investigation  it  was  discovered  that  the  argon  system 
was  much  stiffer  and  a  new  ODE  solver  was  required  for  an  efficient  and  accurate  solution.  The 
Radau5  ODE  integration  solver  was  selected.  This  solver  uses  a  3-stage  5th  order  implicit  Runge- 
Kutta  scheme  with  adaptive  time  stepping.  This  method  allows  for  an  optimum  time  step  to  be 
selected,  which  was  an  improvement  over  the  constant  time  step  that  the  backwards  Euler  method 
used. 


Figure  3:  Example  argon  spectra. 


Another  improvement  made  to  the  CR  model  is  the  ability  to  compute  a  line-by-line  radiation 
spectra  from  the  CR  rates,  as  shown  in  Figure  3.  The  current  algorithm  uses  the  CR  model  to 
generate  the  non-equilibrium  population  of  the  excited  states,  this  in  turn  is  used  to  determine  the 
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1  ELASTIC  COLLISIONS  AND  CR  MODELS 


T  [eV]  T  [eV] 


(a)  Excitation  (2  — >  5)  (b)  Deexcitation  (2  < —  5) 

Figure  4:  Comparison  of  zeroth-order  reaction  rate  with  Monte  Carlo  integration  of  the  full  transfer  integral,  where 
A  =  -p2-.  A  higher  A  indicates  a  larger  disparity  between  the  kinetic  energy  between  the  two  species. 


emission  and  absorption  coefficients.  From  this  data  the  radiation  spectra  is  generated  and  a  plot 
of  intensity  versus  wave  number  is  generated.  This  spectra  data  serves  a  very  useful  validation  tool 
because  the  computationally  generated  spectra  can  be  directly  compared  to  experimental  results. 
Conversely  this  also  provides  experimentalists  with  a  diagnostics  tool  that  can  produce  spectral 
data  for  cases  in  which  the  excited  states  have  not  reached  an  equilibrium.  Currently,  the  model 
assumes  an  optically  thin  plasma,  which  is  only  valid  for  low  density  plasmas,  and  this  is  a  deficiency 
that  should  be  addressed  in  the  future. 

Multi-Fluid  Inelastic  Collisions 

For  this  investigation  the  CR  model  was  extended  to  the  Multi-Fluid  regime.  This  involved  the 
formulation  of  new  exchange  rates  due  to  the  drifting  Maxwellian  approximation.  The  results  of  this 
work  has  resulted  in  two  published  journal  articles;  the  first  covered  excitation/de-excitation [  ],  and 
the  second  includes  ionization  and  recombination  [  ].  The  new  algorithm  is  consistent  with  kinetic 
theory,  has  the  correct  asymptotic  limits,  and  has  microscopic  detailed  balance.  The  excitation/de¬ 
excitation  results  were  compared  to  a  Monte  Carlo  calculation,  shown  in  Figure  4,  and  resulted 
in  a  very  tight  agreement.  A  similar  comparison  was  done  in  the  ionization/recombination  paper 
and  the  Monte  Carlo  integration  had  excellent  agreement  with  the  full  transfer  integral[6].  A  key 
finding  during  this  investigation  was  the  effect  that  the  drift  velocity  has  on  the  kinetics.  Drift 
temperature  is  a  convenient  term  to  express  the  velocity  difference  between  species  and  is  expressed 
as, 

1  msmt 

rp  _  2  ma+mt  1  s 

k 

This  effect  is  shown  in  Figure  5.  This  figure  shows  two  different  results,  one  with  an  initial  drift 
temperature  of  0.01  eV  and  the  other  with  an  initial  drift  temperature  of  10  eV.  For  each  case  the 
simulation  was  run  both  with  and  without  the  exchange  terms  due  to  inelastic  collisions.  For  the 
10  eV  case  the  drift  temperature  is  plotted  in  log-log  scale  beneath  the  semi-log  plot  containing  a 
plot  of  the  other  temperatures.  For  both  cases  the  total  density  is  1020  particles  per  cubic  meter, 
the  initial  ionization  fraction  is  10%  and  the  number  density  for  the  excited  states  is  set  to  a  small 
positive  number,  10-15  particles  per  cubic  meter. 
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2  ADVANCED  PIC  ALGORITHMS 


time  [sec]  time  [sec] 

Figure  5:  Time  evolution  of  the  temperatures  for  two  different  simulations  are  shown.  Each  simulation  was  run  with 
(solid)  and  without  (dashed)  the  exchange  terms  due  to  inelastic  collisions.  The  plot  on  left  starts  with  a  Tw  of  0.01 
eV  and  the  plots  on  the  right  starts  at  Tw  of  10  eV.  Both  simulations  have  a  total  density  of  1020m-'5  and  an  initial 
ionization  fraction  of  10%.  The  plot  on  the  upper  right  shows  the  Boltzmann,  electron,  and  neutral  temperatures 
while  the  plot  on  the  lower  right  shows  the  drift  temperature  in  log-log  scale. 


This  new  algorithm  is  of  particular  relevance  for  laser-produced  plasmas  due  to  the  large  electric 
field  produced  either  directly  from  the  absorption  process  or  via  ponderomotive  forces.  This  E- 
field  creates  an  electron  current  which  could  impose  and/or  sustain  the  relative  drift  between  the 
electrons  and  the  heavy  particles,  which  in  turn  impacts  the  rates.  It  was  also  discovered  that 
the  inelastic  processes  have  an  impact  on  the  overall  thermal  equilibration  process  and  should  be 
included.  These  terms  are  typically  not  included  in  the  standard  collisional  models. 


2  Advanced  PIC  Algorithms 

A  new  method  to  merge  and  split  macro  particles  in  PIC  and  DSMC  codes  was  developed.  This  de¬ 
velopment  started  under  AFOSR  task,  12RZ06COR,  and  was  completed  under  the  current  project. 
The  findings  were  published  in  a  journal  article//.  This  method  was  designed  to  ensure  conser¬ 
vation  and  to  reduce  thermalization  there  by  preserving  the  velocity  distribution  function.  The 
ability  to  control  the  number  of  macro  particles  in  each  cell  is  crucial  to  controlling  the  statistical 
error  both  locally  and  globally.  Two  situations  in  which  this  technique  is  useful  are  when  new 
particles  are  generated  (e.g.  ionization)  or  when  the  density  varies  greatly  throughout  the  domain 
(e.g.  vacuum  plume  expansion). 

As  an  example,  consider  an  ionization  breakdown  simulation  wherein  a  partially  ionized  flow  is 
further  ionized.  A  series  of  x-t  plots  comparing  the  electron  number  density  and  the  computational 
particle  count  both  with  and  without  merging  &  splitting  is  shown  as  labeled  on  the  top  of  Figure  6. 
This  figure  demonstrates  the  schemes  ability  to  control  the  macro  particle  count  while  maintaining 
the  correct  physical  number  density.  Figure  7  shows  the  computation  cost  versus  the  iteration 
count.  When  no  merging  and  splitting  is  used  the  macro  particle  count  grows  at  an  exponential 
rate.  Figures  6  and  7  shows  that  a  significant  time  savings  is  gained  through  the  merge  &  split 
algorithm  without  reducing  the  quality  of  the  results. 
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3  CONSERVATIVE  BGK  INTEGRATOR 
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Figure  6:  Physical  and  computational  electron  density  during  the  ionizing  breakdown  case  with  and  without 
merging  and  splitting  for  a  250  V  potential  from  cathode  (0  mm)  to  anode  (5  mm)  of  each  sub-figure.  The  case 
without  merging  is  labeled  as  control. 


Figure  7:  Computation  cost  during  the  ionizing  breakdown  case  with  and  without  merging  and  splitting  for  a  250 
V  potential  from  cathode  to  anode. 


This  work  has  been  leveraged  by  another  AFOSR  task,  14RQ13COR,  where  this  method  was 
a  critical  component  in  a  fractional  Direct  Simulation  Monte  Carlo  (DSMC)  collision  algorithm. 
The  fractional  particle  collision  algorithm  creates  extra  particles  which  can  make  these  algorithms 
unusable  but  when  combined  with  the  octree  merging  algorithm  the  number  of  particles  and  com¬ 
putational  cost  can  stay  relatively  constant. 


3  Conservative  BGK  Integrator 


One  of  the  difficulties  associated  with  multiscale  simulations  is  the  large  number  of  time  steps 
required  for  completion.  During  the  course  of  a  simulation  errors  can  build  up  and  can  cause 
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numerical  heating  or  cooling.  Traditional  BGK  collisional  operators  exhibit  these  problems  when 
multiple  species  are  allowed  to  collide.  Conservation  becomes  more  of  an  issue  as  the  mass  ratio 
between  the  species  increases,  e.g.  ion/electron  systems.  These  issues  led  to  the  development  of  a 
new  BGK  collision  operator.  This  method  is  fully  conservative  even  when  taking  a  time  step  larger 
than  the  collisional  frequencies.  A  white  paper  containing  the  derivation  and  additional  examples 
are  included  in  Appendix  A. 

To  test  this  algorithm,  the  relaxation  of  a  plasma  consisting  of  hot  electrons  and  cold  ions  is  shown 
in  Figure  8.  In  this  figure  the  same  initial  conditions  are  used  but  the  problem  is  solved  with 
different  time  steps.  Although  the  rate  of  relaxation  is  different,  the  two  solvers  result  in  the 
correct  equilibrium  temperature. 


Figure  8:  Impact  of  time  integration  scheme  for  BGK  equilibration  of  electron-Ar+  plasma. 


This  method  could  benefit  through  the  addition  of  a  velocity  dependent  collision  frequency  or  by 
extending  the  algorithm  to  high  order. 

4  Multiscale  Hybridization  Techniques 

The  last  topic  investigated  was  the  use  hybridization  techniques  to  address  multiscale  systems. 
When  solving  a  physical  system,  there  are  typically  multiple  ways  to  model  the  physical  system 
and  each  model  can  be  solved  using  a  variety  of  numerical  algorithms.  The  techniques  described  in 
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this  section  seek  to  actively  switch  between  different  physical  models;  and/or  numerical  algorithms 
to  efficiently  solve  a  physical  system.  This  new  solver  is  intended  to  be  more  accurate  than  a  coarse 
solution  but  not  as  numerically  expensive  as  a  fine  grained  solution.  A  paper  detailing  this  method 
has  been  written  and  submitted  to  a  journal  for  review. 

For  example,  consider  the  motion  of  a  charged  particle  in  a  magnetic  field.  Figure  9  solves  this 
system  using  four  different  methods,  Gyrokinetics,  leapfrog  PIC,  a  simple  hybrid  scheme,  and  by  a 
modified  time  parallel  method.  The  Gyrokinetics  solver  provides  a  coarse  solution  to  the  particle 
motion;  it  is  computationally  efficient  but  is  expected  to  produce  larger  errors.  The  leapfrog  PIC 
method  requires  a  small  time  step  for  stability  but  produces  a  much  more  accurate  solution  and  is 
used  as  our  fine  solution.  The  hybridization  schemes  represents  two  different  physical  models  that 
governs  the  motion  of  the  particle.  The  simple  hybrid  method  uses  an  error  estimator  to  determine 
which  algorithm  will  be  used  to  solve  the  system  in  the  next  time  step.  The  time  parallel  method 
also  uses  an  error  estimator  to  select  the  solution  method  but  also  provides  a  blending  zone  wherein 
both  solutions  are  solved  1 .  Figure  9  plots  the  relative  error  compared  to  the  computation  time 
required  to  complete  the  simulation.  For  a  given  error  level,  the  plot  towards  the  left  side  represents 
the  most  computationally  efficient  method.  Typically,  a  reduced  time  step  provides  more  accurate 
solutions  at  the  expense  of  increased  computational  time.  The  Gyrokinetics  method  has  in  the 
shortest  computation  time  but  its  error  doesn’t  improve  beyond  a  certain  level.  This  abrupt  end 
to  its  accuracy  is  due  to  the  physical  limitations  of  the  Gyrokinetic  method.  Though  the  leapfrog 
method  converges  to  the  lowest  relative  error,  it  is  ten  to  a  hundred  times  more  computationally 
expensive  than  the  time  parallel  method  in  the  regions  that  the  time  parallel  method  was  used. 
The  simple  hybrid  method  represents  an  improvement  over  the  leapfrog  method  when  the  error  is 
high  but,  just  like  the  Gyrokinetic  method,  fails  to  improve  as  the  computational  effort  increase. 
The  time  parallel  method  represents  a  good  hybridization  of  the  coarse  and  fine  solvers.  For  a  given 
error  it  is  generally  more  accurate  than  the  coarse  solution  and  is  always  more  computationally 
efficient  than  the  fine  solution. 

1Note  that  although  the  time-parallel  method  can  run  in  parallel  all  simulations  were  performed  in  serial. 


Relative  Error 
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Figure  9:  Relative  error  when  computing  the  motion  of  a  changed  particle  in  a  magnetic  mirror 


Moving  forward,  this  technique  will  be  applied  to  other  problems  of  interest  such  as  the  hybridiza¬ 
tion  of  the  Quasi  Steady  State  method  with  the  finite  rate  CR  methods  discussed  earlier. 
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Appendix  A  Conservative  BGK 

A.l  Introduction 

The  Bhatnagar-Gross-Krook  (BGK)  model  is  a  simplified  collisional  model  that  allows  a  particle 
distribution  function  to  relax  towards  a  equilibrium  conditions  at  a  rate  defined  by  a  collisional 
frequency.  Although  not  as  physically  accurate  as  the  full  Boltzmann  collisional  integrals  or  the 
Fokker-Planck  equation,  its  numerical  efficiency  and  simplicity  ensures  its  continued  use  in  multi¬ 
dimensional  kinetic  simulations.  Although  the  scheme  is  conservative  by  nature,  the  majority  of 
the  numerical  implementations  are  unable  to  conserve  the  bulk  properties,  i.e.,  mass,  momentum, 
and  energy,  of  the  gas.  These  numerical  implementations  only  tend  to  preserve  these  quantities 
in  the  limit  of  an  infinitely  small  time  step.  Furthermore,  the  schemes  ability  to  conserve  tend  to 
diminish  as  the  mass  ratio  of  the  gases  increases,  e.g.,  ion-electron. 


A. 2  Derivation 


To  begin  consider  the  standard  multi-component  BGK  operator, 

NS 


£fi 

Dt 


^  ^  Vij  {Fij  fi)  , 


(1) 


where  fi  is  the  velocity  distribution  of  the  ith  species,  Ns  is  the  total  number  of  species  present, 
FtJ  is  the  Maxwellian  of  species  i  with  respect  to  species  j,  and  v1:]  is  the  collisional  rate  of  species  j 
with  species  i.  For  this  investigation  ut]  assumed  to  be  constant  in  velocity  space.  The  Maxwellian 
has  the  form, 
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where  n*  and  r/ij  are  the  number  density  and  the  mass  of  species  i  respectively,  Kb  is  the  Boltzmann 
constant  and  Af  is  the  number  of  velocity  dimensions.  Tij  and  Ujj  are  the  characteristic  temperature 
and  bulk  velocity  in  A I  dimensions.  The  calculation  of  m,  T \j  and  u ij  will  be  detailed  later.  We 
can  rewrite  Eq.  (1)  as, 
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For  simplicity  let  A,  =  "ij  and  Kf>  =  Jfj 

constant  in  time  and  defined  at  tn  then  Eq.  (1)  has  an  analytical  solution, 


y .  If  we  assume  that  Uij  and  F. \j  are 
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where  ff  is  the  distribution  of  fi  at  time  tn  and  A t  =  t  —  tn.  For  simplicity  let 
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As  a  side  note  Eq.  (1)  may  also  be  solved  through  an  Euler  method 

m  =  f? +4](^f 

where 

0  implicit  Euler 

E  = 
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At 


1  +  AtK[i]  (1  -  E) 


1/2  Crank-Nicolson 
1  explicit  Euler 


We  return  our  attention  to  the  bulk  parameters,  nt,  Tjj ,  and  Uy.  For  the  case  when  j  =  i  the  bulk 
parameters  are  determined  by  taking  the  moments  of  the  distribution  function  /j  at  time  tn, 

ni  =  [  f?dv, 


Uu  =  Uj  =  —  /  v/f'dv, 
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The  values  of  Rj  and  Uj  f  when  j  ^  i  are  found  by  enforcing  momentum  and  energy  conservation. 
Momentum  is  conserved  by  enforcing, 
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J  V 


Enforcing  this  over  a  single  time  step  yields, 
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or  substituting  Eq.  (5)  into  Eq.  (8) 
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and  then  integrating  over  velocity 
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Eq.  (10)  represents  one  equation  and  (Ns  —  1  )Ng  unknowns,  so  additional  constraints  are  required. 
We  enforce  that  as  time  advances  all  species  approach  the  same  velocity,  i.e.,  u,j  =  u,  for  j  /  i. 
With  this  added  constraint  Eq.  (10)  can  be  solved  for  u, 
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Rj  is  determined  by  enforcing  energy  conservation  over  a  time  step, 
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By  following  the  same  steps  as  in  the  momentum  conservation  and  by  assuming  that  all  species 
will  equilibrate  towards  the  same  temperature,  Ttj  =  T  for  j  /  i,  yields, 
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Substituting  Eq.  (11)  and  Eq.  (13)  into  Eq.  (2)  yields 
Fij(y,T,  u)  =  m 
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It  should  be  noted  that  applying  conservation  of  mass  does  not  provide  any  useful  relationships 
because  the  only  terms  present  are  the  mass  and  number  density  of  each  species.  This  is  easily 
shown  by  taking  the  zeroth  moment, 

NS  ,  NS 


X>  f  fi(tn  +  A t)dv  =  y,  rrij  f  fi(tn)dv 
■  J  v  J  V 


NS 


NS 


(15) 


'Y^mini  =  Y.mirn. 


This  derivation  was  done  for  an  arbitrary  number  of  species,  but  in  practice  it  is  much  more  common 
to  have  only  one  or  two  specie  present.  As  such  the  following  equations  are  a  simplification  of  the 
above  assuming  that  only  two  species,  a  and  b  are  present. 
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Faa.  Fab,  Fbb ,  and  Fba  are  easily  computed  from  Eq.  (2)  and  Eq.  (14). 

It  is  worthwhile  to  draw  attention  to  several  facts  about  this  new  derivation, 


•  The  selection  of  the  collisional  rates  are  independent  of  each  other,  and  does  not  affect  con¬ 
servation. 

•  No  special  treatment  is  required  if  both  ions  and  electrons  are  present  in  the  plasma  to  ensure 
conservation. 

•  Ensuring  that  mass,  momentum,  and  energy  are  equal  at  tn  and  tn  +  At,  is  a  weaker  require¬ 
ment  than  enforcing  jjj  fv  f(t)v°'1,2dv  =  0,  but  the  later  method  does  not  have  an  analytical 
solution. 

•  The  only  difference  between  the  analytical  solution,  implicit  Euler,  explicit  Euler,  and  Crank- 
Nicolson  methods  is  in  the  initial  evaluation  of  Aa  and  Ae- 

•  Tj  and  are  also  functions  of  time  and  can  be  evaluated  at  tn  +  At  with  out  numerically 
integrating  over  velocity  space. 
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This  section  demonstrates  the  schemes  ability  to  conserve  mass,  momentum,  and  energy  by  applying 
the  new  scheme  to  two  different  test  cases.  The  first  test  consists  relaxation  of  both  temperature 
and  bulk  velocity  for  Argon  and  Xenon  ions  in  three-dimensional  velocity  space.  The  second  case 
shows  temperature  relaxation  between  Argon  ions  and  electrons;  this  test  case  is  performed  with 
one  velocity  dimension.  These  test  cases  will  show  that  mass,  momentum,  and  energy  are  conserved 
to  machine  precision  and  that  the  entropy  of  the  system  increases  in  time.  The  results  will  also 
show  that  the  temperature  and  bulk  velocity  relaxes  towards  the  correct  equilibrium  conditions. 

The  collisional  frequencies  used  in  this  investigation  were  obtained  from  the  NRL  plasma  formulary 
with  a  small  change  for  use  with  the  inks  system  and  to  enforce  a  slower  relaxation  of  the  ions, 
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where  T  is  equal  to  the  temperature  in  eV,  mp  is  the  mass  of  a  proton,  and  z  is  the  unit  charge  of 
the  ion.  uab  is  for  ion-ion  collisions  where  a  ^  b.  The  Coulomb  logarithm  for  all  cases  is  equal  to, 
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where  So  is  the  permittivity  of  free  space  and  qe  is  the  charge  of  an  electron. 


A. 3.1  Argon-Xenon  Thermal  and  Momentum  Equilibrium 

For  this  test  case  we  consider  the  following  initial  conditions  for  Argon:  n  =  2.5  x  1019  (1/m3), 
T  =  1000  (Kelvin),  u  =  (1000,0,0)  (m/s),  and  z  =  1  (electron  charge).  And  Xenon:  n  = 
2.5  x  1019  (1/m3),  T  =  10000  (Kelvin),  u  =  (—1000,0,0)  (m/s),  and  z  =  1  (electron  charge). 
Velocity  space  is  bounded  at  ±2.5  x  104  (nr/s)  in  all  directions  and  the  mesh  spacing  is  100  m/s 
in  all  directions.  The  simulation  ran  for  5  micro-seconds  with  a  time  step  of  10  nano-seconds.  The 
velocity  bounds  and  spacing  were  selected  such  that  the  correct  bulk  properties,  number  density, 
velocity  and  temperature,  were  obtained  at  the  start  of  the  simulation.  The  equilibrium  conditions 
are  determined  by  enforcing  momentum  and  energy  conservation, 
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where  the  superscript  0  represents  the  initial  conditions  and  the  subscript  eq  represents  the  equi¬ 
librium  conditions.  With  these  initial  conditions  the  equilibrium  velocity  and  temperature  are, 
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Figure  A-l:  Velocity  and  temperature  time  history  for  Argon- Xenon  equilibration. 


-533.421320034  m/s  and  7955.8852308  Kelvin  respectively.  Figure  A-l  shows  that  the  correct 
equilibrium  conditions  are  achieved,  while  the  results  shown  in  Fig.  A- 2  show  that  mass, 
momentum,  and  energy  are  conserved  to  machine  precision  and  that  entropy  increases. 


A. 3. 2  Argon  Electron  Thermal  Equilibrium 

For  this  test  case  we  consider  the  following  initial  conditions  for  Argon:  n  =  2.5  x  1019  (1/m3), 
T  =  1000  (Kelvin),  u  =  0  (m/s),  and  z  =  1  (electron  charge).  Velocity  space  for  Argon  is  bounded 
at  ±5  x  104  m/s  with  a  spacing  of  250  m/s.  The  electrons  have  the  following  initial  conditions, 
n  =  2.5  x  1019  (1/m3),  T  =  10000  (Kelvin),  and  u  =  0  (m/s).  Velocity  space  for  the  electrons 
is  bounded  at  ±1  x  107  m/s  with  a  spacing  of  50000  m/s.  The  simulation  was  performed  using 
several  time  steps  ranging  from  10_8to  10_11and  vas  allowed  to  run  f  or  3x  10~4seconds.  Atypical 
simulation  isshown  inRgsA-3,  and  A-4,  and  used  a  time  step  of  10“10  seconds.  In  Fig.  A-3  the  time 
history  of  velocity  and  temperature  for  the  Argon  electron  equilibration  is  shown;  as  expected  the 
velocity  of  the  two  species  stays  near  zero  and  the  temperatures  of  the  two  gases  approaches  the 
correct  equilibrium  temperature  of  5500  Kelvin.  In  Fig.  A-4  the  time  history  of  the  conservation  bulk 
properties  as  well  as  the  entropy  are  shown.  The  most  likely  cause  for  this  is  an  accumulation  of 
round  off  error,  which  is  supported  by  the  fact  that  the  error  per  time  step  is  around  machine 
precision. 

Even  though  the  bulk  properties  are  preserved  the  accuracy  of  the  temperature  and  velocity  profiles 
are  affected  by  the  time  step.  Figure  A-5  shows  the  time  history  of  the  error  in  the  ion  temperature  as 
well  as  the  convergence  rate  of  temperature.  This  plot  shows  the  temperature  convergence  at  a  first 
order  rate.  This  is  not  unexpected  since  it  is  assumed  that  the  collisional  frequency  and  the 
Maxwellian  are  constant  during  a  time  step,  effectively  a  first  order  Taylor  series  expansion  in  time. 
It  should  be  possible  to  increase  the  order  of  accuracy  by  expressing  the  collisional  frequency  and  the 
Maxwellian  as  a  higher  order  Taylor  series  in  time.  Also  of  note  is  that  no  analytical  solution  exists 
for  this  problem  and  the  simulation  was  not  run  to  convergence,  thus  the  reference  solution  was 
chosen  at  a  time  step  of  10-11  seconds. 
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Figure  A-2:  Time  history  plot  of  the  conservation  properties  and  entropy  for  Argon- Xenon 
equilibration. 
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Figure  A-3:  Velocity  and  temperature  time  history  for  Argon-electron  equilibration. 
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Figure  A-4:  Time  history  plot  of  the  conservation  properties  and  entropy  for  Argon-electron 
equilibration.  For  conservation  plots  the  top  is  cumulative  and  bottom  is  per  time  step. 


(a)  Relative  error  in  ion  temperature 
Figure  A-5:  Velocity  and  temperature  time  history  for  Argon-electron  equilibration. 
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